## Callin Switzer
## 17 Nov 2017
## Multilevel model to visualize bees'
## behavior on the artificial pollen system

#install packages
ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c("ggplot2", "reshape2", 'lme4', 'sjPlot', "multcomp", "plyr", "effects")
ipak(packages)
## Loading required package: ggplot2
## Loading required package: reshape2
## Loading required package: lme4
## Loading required package: Matrix
## Loading required package: sjPlot
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.10
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source or restore original 'Matrix' package
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
## Loading required package: plyr
## Loading required package: effects
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
##  ggplot2 reshape2     lme4   sjPlot multcomp     plyr  effects 
##     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE
# set ggplot theme
theme_set(theme_classic() + theme(axis.text=element_text(colour="black")))

# define data and figure directories
dataDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehData"
figDir <- "/Users/cswitzer/Dropbox/SonicationBehavior/SonBehFigs"

print(paste("last run ", Sys.time()))
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2017c.1.0/
## zoneinfo/America/Los_Angeles'
## [1] "last run  2017-12-08 22:58:40"
print(R.version)
##                _                           
## platform       x86_64-apple-darwin15.6.0   
## arch           x86_64                      
## os             darwin15.6.0                
## system         x86_64, darwin15.6.0        
## status                                     
## major          3                           
## minor          4.2                         
## year           2017                        
## month          09                          
## day            28                          
## svn rev        73368                       
## language       R                           
## version.string R version 3.4.2 (2017-09-28)
## nickname       Short Summer

Read in data and double check it

sl <- read.csv(file.path(dataDir, '01_CombinedTrials_cleaned.csv'))

# make sure hive is a factor
sl$hive = as.factor(sl$hive)

head(sl)
##   beeCol hive meanFreq IT_imputed index freq     amp
## 1   blue    4 395.2778       3.61     1  400 0.06328
## 2   blue    4 395.2778       3.61     2  360 0.34299
## 3   blue    4 395.2778       3.61     3  430 0.44099
## 4   blue    4 395.2778       3.61     4  420 0.16841
## 5   blue    4 395.2778       3.61     5  420 0.15284
## 6   blue    4 395.2778       3.61     6  410 0.18302
##                    datetime rewNum rewTF lowRewAmp highrewAmp
## 1  2016_11_22__08_23_47_721      1     T         0          5
## 2  2016_11_22__08_23_49_677      2     T         0          5
## 3  2016_11_22__08_23_50_139      3     T         0          5
## 4  2016_11_22__08_24_10_938      4     T         0          5
## 5  2016_11_22__08_24_18_013      5     T         0          5
## 6  2016_11_22__08_24_18_464      6     T         0          5
##                         BeeNumCol
## 1 BeeBlue_22Nov2016_Hive4_initial
## 2 BeeBlue_22Nov2016_Hive4_initial
## 3 BeeBlue_22Nov2016_Hive4_initial
## 4 BeeBlue_22Nov2016_Hive4_initial
## 5 BeeBlue_22Nov2016_Hive4_initial
## 6 BeeBlue_22Nov2016_Hive4_initial
##                                     accFile trialNum
## 1 2016_11_22__08_23_47_721_220_450_test.txt        1
## 2 2016_11_22__08_23_49_677_220_450_test.txt        1
## 3 2016_11_22__08_23_50_139_220_450_test.txt        1
## 4 2016_11_22__08_24_10_938_220_450_test.txt        1
## 5 2016_11_22__08_24_18_013_220_450_test.txt        1
## 6 2016_11_22__08_24_18_464_220_450_test.txt        1
##                 datetime_str lowFrq highFrq   IT  trt   amp_acc
## 1 2016-11-22 08:23:47.721000    220     450 3.61 full  6.222222
## 2 2016-11-22 08:23:49.677000    220     450 3.61 full 33.725664
## 3 2016-11-22 08:23:50.139000    220     450 3.61 full 43.361849
## 4 2016-11-22 08:24:10.938000    220     450 3.61 full 16.559489
## 5 2016-11-22 08:24:18.013000    220     450 3.61 full 15.028515
## 6 2016-11-22 08:24:18.464000    220     450 3.61 full 17.996067
dim(sl) # should be 24303 rows
## [1] 24303    21
# hive 5 is most common
table(sl$hive, useNA = 'always')
## 
##     3     4     5  <NA> 
##   513  1047 22743     0
## make sure all bee colors are lowercase
sl$beeCol <- tolower(sl$beeCol)

# make sure there are values lower than 220 and higher than 450 
# (the cutoff for buzzes used in the experiment)
hist(sl$freq, breaks = seq(215, 450, by = 10))

nrow(sl[sl$freq < 220 | sl$freq > 450,]) # should have 0 rows
## [1] 0
# look at treatments
xtabs(~sl$beeCol+ trt, data = sl )
##                   trt
## sl$beeCol          full high  low
##   blue               36    0    0
##   gold               54    0    0
##   goldred             3    0    0
##   green              13    0    0
##   lime               28    0    0
##   limeblue          530    0    0
##   limegold           54    0    0
##   limegreen          50    0    0
##   limeorange         84    0    0
##   limepink           51    0    0
##   limepurple        109    0  243
##   limepurpleyellow   58    0    0
##   limered            50    0 3424
##   limesilver          3    0    0
##   limewhite         101    0    0
##   limeyellow         52    0    0
##   orange             85  416    0
##   orangeblue         50    0    0
##   orangegreen        50    0    0
##   orangepink         50    0    0
##   orangepurple       50    0    0
##   purple              9    0    0
##   redblue           673    0    0
##   redgreen          530    0    0
##   redpink            56    0  621
##   redpurple          55  163    0
##   silver             26    0    0
##   white             283    0    0
##   whiteblue          56    0 2730
##   whitegold          48    0 1177
##   whitegreen         39    0    0
##   whiteorange        54    0 1068
##   whitepink         135 1049    0
##   whitepurple        55    0    0
##   whitered           59 1285    0
##   whiteyellow       157    0    0
##   yellowblue         54 1533    0
##   yellowgreen       532    0    0
##   yelloworange       54 2059    0
##   yellowpink         58    0 2251
##   yellowpurple       54    0 1189
##   yellowred         547    0    0
# find percentage reward by treatment
mean(grepl("[tT]", as.character(sl$rewTF))) # overall mean
## [1] 0.4155454
# percentage that were rewarded by treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, mean)
##      full      high       low 
## 1.0000000 0.3535742 0.2128631
# total number of trials for each treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, length)
##  full  high   low 
##  5095  6505 12703
# total number of rewards per treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, FUN = function(x) sum(x))
## full high  low 
## 5095 2300 2704
# total number of trials that were unrewarded per treatment
tapply((grepl("[tT]", as.character(sl$rewTF))), INDEX= sl$trt, FUN = function(x) sum(!x))
## full high  low 
##    0 4205 9999

Visualizations and modeling

# calculate trial averages and plot
sl$colNum = paste(sl$beeCol, sl$hive, sep = "_")
sl$trt <- as.character(sl$trt)
sl$trt[sl$trt == "full" & sl$trialNum >1 ] <- "full_2"


aggdata <- aggregate(sl$freq, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(aggdata)[colnames(aggdata) == "x"] = "freq"

aggdata_sd <- aggregate(sl$freq, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=sd, na.rm=TRUE)
colnames(aggdata_sd)[colnames(aggdata_sd) == "x"] = "freq_sd"

aggdata = merge(aggdata, aggdata_sd)
#aggdata$trt = sapply(1:nrow(aggdata), FUN = function(x){sl[sl$colNum == aggdata$colNum[x] & sl$trialNum == aggdata$trialNum[x], "trt"][1]})
aggdata = aggdata[order(aggdata$colNum, aggdata$trialNum, decreasing = FALSE), ]
agg_sm = aggdata[aggdata$trialNum <= 2, ]
rownames(agg_sm) = 1:nrow(agg_sm)
agg_sm
##                colNum trialNum    trt     freq  freq_sd
## 1              blue_4        1   full 395.2778 39.38717
## 2              gold_3        1   full 348.3333 37.15089
## 3           goldred_4        1   full 406.6667 32.14550
## 4             green_4        1   full 318.8889 57.32461
## 5             green_4        2 full_2 285.0000 59.16080
## 6              lime_5        1   full 374.6429 36.36055
## 7          limeblue_5        1   full 302.4000 34.73280
## 8          limeblue_5        2 full_2 299.8148 39.02059
## 9          limegold_5        1   full 344.6296 44.07168
## 10        limegreen_5        1   full 370.8000 29.40498
## 11       limeorange_5        1   full 285.7576 35.26953
## 12       limeorange_5        2 full_2 344.5098 26.70683
## 13         limepink_5        1   full 308.4314 43.23760
## 14       limepurple_5        1   full 352.2642 38.66232
## 15       limepurple_5        2    low 344.2798 32.55948
## 16 limepurpleyellow_5        1   full 336.7241 50.76219
## 17          limered_5        1   full 354.6000 37.91559
## 18          limered_5        2    low 365.3381 36.20344
## 19       limesilver_5        1   full 306.6667 11.54701
## 20        limewhite_5        1   full 292.0000 44.90352
## 21        limewhite_5        2 full_2 327.0588 35.51305
## 22       limeyellow_5        1   full 339.4231 39.37818
## 23           orange_3        1   full 388.5294 38.22837
## 24           orange_5        1   full 345.0980 33.60789
## 25       orangeblue_5        1   full 301.4000 45.17585
## 26      orangegreen_5        1   full 296.2000 40.45103
## 27       orangepink_5        1   full 286.6000 27.15113
## 28     orangepurple_5        1   full 329.0000 25.65469
## 29           purple_3        1   full 382.2222 21.08185
## 30          redblue_4        1   full 342.2449 40.53113
## 31          redblue_4        2 full_2 335.6667 51.89189
## 32         redgreen_5        1   full 334.7059 34.19666
## 33         redgreen_5        2 full_2 314.5283 31.53545
## 34          redpink_5        1   full 288.3929 41.06654
## 35          redpink_5        2    low 275.3448 34.24450
## 36        redpurple_5        1   full 311.4545 31.64726
## 37        redpurple_5        2   high 336.4444 40.75995
## 38           silver_5        1   full 305.3846 38.39070
## 39            white_4        1   full 343.1250 47.49720
## 40            white_4        2 full_2 342.2388 38.95689
## 41        whiteblue_5        1   full 318.2143 27.17667
## 42        whiteblue_5        2    low 351.6242 29.13402
## 43        whitegold_5        1   full 315.2083 34.20772
## 44        whitegold_5        2    low 307.2727 45.76009
## 45       whitegreen_4        1   full 300.7692 33.43372
## 46      whiteorange_5        1   full 324.4444 51.71280
## 47      whiteorange_5        2    low 355.2308 44.74234
## 48        whitepink_5        1   full 336.8750 42.68285
## 49        whitepink_5        2   high 323.9161 28.55770
## 50      whitepurple_5        1   full 328.1818 43.50781
## 51         whitered_5        1   full 288.8136 36.05956
## 52         whitered_5        2   high 304.6667 42.10035
## 53      whiteyellow_5        1   full 303.4545 39.07344
## 54      whiteyellow_5        2 full_2 312.5532 36.56252
## 55       yellowblue_5        1   full 313.7037 44.86021
## 56       yellowblue_5        2   high 327.9688 49.76099
## 57      yellowgreen_5        1   full 298.0392 40.64577
## 58      yellowgreen_5        2 full_2 300.6000 27.58364
## 59     yelloworange_5        1   full 268.1481 36.18996
## 60     yelloworange_5        2   high 313.9205 39.62646
## 61       yellowpink_5        1   full 323.6207 39.98676
## 62       yellowpink_5        2    low 335.7042 41.00507
## 63     yellowpurple_5        1   full 354.4444 36.37678
## 64     yellowpurple_5        2    low 373.0058 34.17237
## 65        yellowred_5        1   full 320.1754 41.85396
## 66        yellowred_5        2 full_2 333.0357 29.96481
agg_sm[duplicated(data.frame(agg_sm$colNum, agg_sm$trialNum)), ]
## [1] colNum   trialNum trt      freq     freq_sd 
## <0 rows> (or 0-length row.names)
ggplot(agg_sm, aes(x = trt, y = freq, fill = trialNum > 1)) +
  geom_boxplot(alpha = 0.2) + 
  geom_point(aes(color = trialNum>1)) + 
  geom_line(aes(group = colNum))

diffdf <- sapply(unique(agg_sm$colNum), FUN = function(x){
  tmp = agg_sm[agg_sm$colNum == x, ]
  if(nrow(tmp) <= 1)
    diff = NA
  else
    diff = tmp$freq[tmp$trialNum == 2] - tmp$freq[tmp$trialNum == 1]
  return(diff)
})

trtDF = sapply(unique(agg_sm$colNum), FUN = function(x){
  tmp = agg_sm[agg_sm$colNum == x, ]
  ttrs = paste(tmp$trt[tmp$trialNum == 1], tmp$trt[tmp$trialNum == 2], sep = "_")
  return(ttrs)
})

buzzdiffs = data.frame(trtDF, diffdf)

ggplot(buzzdiffs, aes(x = trtDF, y= diffdf)) + 
  geom_boxplot() + 
  geom_point()
## Warning: Removed 20 rows containing non-finite values (stat_boxplot).
## Warning: Removed 20 rows containing missing values (geom_point).

agg2 = aggregate(sl$freq, by=list(colNum = sl$colNum, fullTrt = sl$trt == "full", trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(agg2)[colnames(agg2) == "x"] = "freq"
agg2$trt = as.character(agg2$trt)

diffdf <- t(as.data.frame(t(sapply(unique(agg2$colNum), FUN = function(x){
  tmp = agg2[agg2$colNum == x, ]
  if(nrow(tmp) <= 1)
    return(NA)
  if (length(unique(tmp$trt)) > 2){
    tmp = tmp[tmp$trt != "full_2", ]
  }
  diff =  tmp$freq[!tmp$fullTrt] - tmp$freq[tmp$fullTrt]
  return(diff)

}))))

length(diffdf)
## [1] 43
trtDF = sapply(unique(agg2$colNum), FUN = function(x){
  tmp = agg2[agg2$colNum == x, ]
  if (length(unique(tmp$trt)) > 2){
    tmp = tmp[tmp$trt != "full_2", ]
  }
  ttrs = paste(tmp$trt[tmp$fullTrt], tmp$trt[!tmp$fullTrt], sep = "_")
  return(ttrs)
})

length(trtDF)
## [1] 43
buzzdiffs = data.frame(trtDF, diffdf)
tapply(buzzdiffs$diffdf, INDEX = buzzdiffs$trtDF, mean)
##       full_ full_full_2   full_high    full_low 
##          NA    4.851689   13.209001   14.307123
ggplot(droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]), aes(x = trtDF, y= as.numeric(diffdf))) + 
  geom_boxplot() + 
  geom_point() + 
  labs(y = "diff in frequency -- averages")

m1 = lm(as.numeric(diffdf) ~ trtDF - 1, data = droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]))
summary(m1)
## 
## Call:
## lm(formula = as.numeric(diffdf) ~ trtDF - 1, data = droplevels(buzzdiffs[buzzdiffs$trtDF != 
##     "full_", ]))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -38.741 -17.992  -3.556  15.767  53.901 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## trtDFfull_full_2    4.852      7.488   0.648    0.524
## trtDFfull_high     13.209      9.667   1.366    0.186
## trtDFfull_low      14.307      8.372   1.709    0.102
## 
## Residual standard error: 23.68 on 21 degrees of freedom
## Multiple R-squared:  0.1987, Adjusted R-squared:  0.08422 
## F-statistic: 1.736 on 3 and 21 DF,  p-value: 0.1904
sl$trt = relevel(factor(sl$trt), ref = "full")


sl$trt2 = mapvalues(sl$trt, from = c("full", "full_2", "high", "low"), 
                    to = c("full", "full", "high", "low"))


# summary for paper
# fit a varying slope and intercept for colNum (bee ID), and allow the slope of the trialNum
# variable to vary by colNum (beeID)
m2 = lmer(freq ~ trt2 + hive + trialNum + IT_imputed + (1+trialNum|colNum), data = sl, REML = FALSE)

# this model estimates a global intercept
# random effect intercept for colNum (beeID)
# a single global estimate for trialNum
# the effect of trialNum within each level of colNum
# the correlation between intercept of trialNum across levels of colNum

m2_1 = lmer(freq ~ trt2* IT_imputed + hive + trialNum  + (1+trialNum|colNum), data = sl, REML = FALSE)

# compare models with BIC
BIC(m2, m2_1) # m2 (no interaction) is better
##      df      BIC
## m2   11 246086.5
## m2_1 13 246090.2
anova(m2, m2_1) # Note: this disagrees with BIC
## Data: sl
## Models:
## m2: freq ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## m2:     colNum)
## m2_1: freq ~ trt2 * IT_imputed + hive + trialNum + (1 + trialNum | 
## m2_1:     colNum)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m2   11 245997 246087 -122988   245975                             
## m2_1 13 245985 246090 -122979   245959 16.547      2  0.0002552 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 = update(m2, .~. - hive)
BIC(m2, m3) # m3 is better (without hive)
##    df      BIC
## m2 11 246086.5
## m3  9 246073.9
anova(m2, m3) # Note: disagrees with BIC
## Data: sl
## Models:
## m3: freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | colNum)
## m2: freq ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## m2:     colNum)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m3  9 246001 246074 -122991   245983                           
## m2 11 245997 246087 -122988   245975 7.5417      2    0.02303 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit best model with REML = TRUE
m3 <- update(m3, .~., REML = TRUE)
summary(m3) # summary for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | colNum)
##    Data: sl
## 
## REML criterion at convergence: 245962
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0242 -0.5161  0.1663  0.6809  3.9135 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  colNum   (Intercept) 1072     32.74         
##           trialNum     112     10.58    -0.73
##  Residual             1439     37.93         
## Number of obs: 24303, groups:  colNum, 43
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  479.138     44.361  10.801
## trt2high      11.161      2.294   4.865
## trt2low       17.718      1.916   9.248
## trialNum       1.448      2.133   0.679
## IT_imputed   -37.291     10.773  -3.462
## 
## Correlation of Fixed Effects:
##            (Intr) trt2hg trt2lw trilNm
## trt2high    0.069                     
## trt2low    -0.005  0.009              
## trialNum    0.007 -0.042 -0.037       
## IT_imputed -0.993 -0.072  0.000 -0.093
plot(m3)

qqnorm(ranef(m3)$colNum[[1]])
qqline(ranef(m3)$colNum[[1]])

sjp.lmer(m3) # plot random effects to find any outliers
## Plotting random effects...

sjp.lmer(m3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).

# post-hoc tests -- pvals for paper
summary(glht(m3, linfct = mcp(trt2 = "Tukey")), test = adjusted("none"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = freq ~ trt2 + trialNum + IT_imputed + (1 + trialNum | 
##     colNum), data = sl, REML = TRUE)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## high - full == 0   11.161      2.294   4.865 1.14e-06 ***
## low - full == 0    17.718      1.916   9.248  < 2e-16 ***
## low - high == 0     6.557      2.975   2.204   0.0275 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)

Generate CI’s

# set number of bootstrap samples
nbootSims = 1000

ITmean = mean(tapply(sl$IT_imputed, INDEX = sl$beeCol, FUN = function(x) x[1] ))

# don't need hive, because that's not in the model we chose (above)
pframe <- data.frame(trt2 = levels(droplevels(sl$trt2)), IT_imputed = ITmean,  colNum = 99999, trialNum = 2)
pframe$freq <- 0
pp <- predict(m3, newdata = pframe, re.form=NA, type = 'response') # re.form sets all random effects to 0


### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(m3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))
pframe$blo<-bb2_se[1,]
pframe$bhi<-bb2_se[2,]
pframe$predMean <- pp
pframe <- pframe[, c('trt2', "blo", "bhi", "predMean")]
pframe
##   trt2      blo      bhi predMean
## 1 full 322.2094 336.6227 329.4150
## 2 high 332.2079 349.1315 340.5761
## 3  low 339.6179 354.8610 347.1330

Make plots for paper

#plot 95% confidence intervals
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
# rename labels for plot
pframe$trt3 = mapvalues(pframe$trt2, from = c("full", "high", "low"), 
                        to = c("Full range\n(220 - 450 Hz)", "High range\n(340 - 390 Hz)", "Low range\n(220 - 330 Hz)"))

pframe
##   trt2      blo      bhi predMean                       trt3
## 1 full 322.2094 336.6227 329.4150 Full range\n(220 - 450 Hz)
## 2 high 332.2079 349.1315 340.5761 High range\n(340 - 390 Hz)
## 3  low 339.6179 354.8610 347.1330  Low range\n(220 - 330 Hz)
g0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
     geom_point()+ 
     labs(y = "Sonication frequency (Hz)", x = "Frequency range for reward") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
  annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "c"),
                color="black") 
g0

ggsave(plot = g0, filename = file.path(figDir, "SonicationFreqPredsAndCI_unadjusted.pdf"), width = 5, height = 4)



g0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
     geom_point()+ 
     labs(y = "Sonication frequency (Hz)", x = "Frequency range for reward") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
  annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "b"),
                color="black") 
g0

ggsave(plot = g0, filename = file.path(figDir, "SonicationFreqPredsAndCI_adjustedPvals.pdf"), width = 5, height = 4)



# use the "effects" library to get CI's
# this basically gives the same result as bootstrapping
# it's a good double check
ee <- as.data.frame(Effect(c("trt2"), 
                           fixed.predictors =list(given.values = c(trialNum = 2, 
                                                                   IT_imputed = ITmean)),  
                           m3) )


# compare two methods -- very similar
ee
##   trt2      fit       se    lower    upper
## 1 full 329.4150 3.798940 321.9689 336.8612
## 2 high 340.5761 4.280727 332.1856 348.9666
## 3  low 347.1330 4.081495 339.1330 355.1330
pframe
##   trt2      blo      bhi predMean                       trt3
## 1 full 322.2094 336.6227 329.4150 Full range\n(220 - 450 Hz)
## 2 high 332.2079 349.1315 340.5761 High range\n(340 - 390 Hz)
## 3  low 339.6179 354.8610 347.1330  Low range\n(220 - 330 Hz)
g1 <- ggplot(ee, aes(x=trt2, y=fit))+
     geom_point()+ 
     labs(y = "Sonication Frequency (Hz)", x = "Treatment") + 
     geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1)+
     theme_classic() + 
    annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 360, label=c("a", "b", "c"),
                color="black") 
g1

Analysis for amplitude

# refref: here try log amp_acc and REML = FALSE with BIC

# interaction model
# note log-transformation to make model fit assumptions better
maa0 = lmer(log(amp_acc) ~ trt2* IT_imputed + hive + trialNum  + (1+trialNum|colNum), data = sl, REML = FALSE)

# main effect model
maa1 = lmer(log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1+trialNum|colNum), data = sl, REML = FALSE)
BIC(maa0, maa1) # use no interaction (keep maa1)
##      df      BIC
## maa0 13 44377.70
## maa1 11 44360.43
anova(maa0, maa1) # note that this agrees with likelihood ratio test
## Data: sl
## Models:
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## maa1:     colNum)
## maa0: log(amp_acc) ~ trt2 * IT_imputed + hive + trialNum + (1 + trialNum | 
## maa0:     colNum)
##      Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## maa1 11 44271 44360 -22125    44249                         
## maa0 13 44272 44378 -22123    44246 2.9209      2     0.2321
maa2 = update(maa1, .~. - trt2)
BIC(maa1, maa2) # keep treatment (maa1)
##      df      BIC
## maa1 11 44360.43
## maa2  9 44375.90
anova(maa1, maa2) # agrees with LRT; p-val for paper, if needed
## Data: sl
## Models:
## maa2: log(amp_acc) ~ hive + trialNum + IT_imputed + (1 + trialNum | 
## maa2:     colNum)
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## maa1:     colNum)
##      Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
## maa2  9 44303 44376 -22142    44285                             
## maa1 11 44271 44360 -22125    44249 35.665      2  1.801e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
maa3 <- update(maa1, .~. - hive)
BIC(maa1, maa3) # remove hive (maa3)
##      df      BIC
## maa1 11 44360.43
## maa3  9 44345.85
anova(maa1, maa3) # agrees with LRT
## Data: sl
## Models:
## maa3: log(amp_acc) ~ trt2 + trialNum + IT_imputed + (1 + trialNum | 
## maa3:     colNum)
## maa1: log(amp_acc) ~ trt2 + hive + trialNum + IT_imputed + (1 + trialNum | 
## maa1:     colNum)
##      Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)  
## maa3  9 44273 44346 -22128    44255                           
## maa1 11 44271 44360 -22125    44249 5.6143      2    0.06038 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit with REML = TRUE for paper
maa3 <- update(maa3, .~., REML = TRUE)
summary(maa3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ trt2 + trialNum + IT_imputed + (1 + trialNum |  
##     colNum)
##    Data: sl
## 
## REML criterion at convergence: 44279.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8377 -0.5961  0.1156  0.7367  2.9119 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  colNum   (Intercept) 0.07079  0.26607       
##           trialNum    0.00203  0.04506  -0.57
##  Residual             0.35872  0.59893       
## Number of obs: 24303, groups:  colNum, 43
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 2.398060   0.455225   5.268
## trt2high    0.194526   0.034446   5.647
## trt2low     0.067759   0.029326   2.311
## trialNum    0.005513   0.010107   0.545
## IT_imputed  0.324122   0.110268   2.939
## 
## Correlation of Fixed Effects:
##            (Intr) trt2hg trt2lw trilNm
## trt2high    0.107                     
## trt2low    -0.009  0.029              
## trialNum    0.036 -0.089 -0.083       
## IT_imputed -0.995 -0.114  0.000 -0.087
# diagnostics
plot(maa3)

qqnorm(ranef(maa1)$colNum[[1]])
qqline(ranef(maa1)$colNum[[1]])

sjp.lmer(maa3) # plot random effects to find any outliers
## Plotting random effects...

sjp.lmer(maa3, type = 'fe', sort = TRUE, p.kr = FALSE) # plot fixed effects
## Computing p-values via Wald-statistics approximation (treating t as Wald z).

# post-hoc tests -- pvals for paper
summary(glht(maa3, linfct = mcp(trt2 = "Tukey")), test = adjusted("none"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = log(amp_acc) ~ trt2 + trialNum + IT_imputed + 
##     (1 + trialNum | colNum), data = sl, REML = TRUE)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## high - full == 0  0.19453    0.03445   5.647 1.63e-08 ***
## low - full == 0   0.06776    0.02933   2.311  0.02086 *  
## low - high == 0  -0.12677    0.04459  -2.843  0.00447 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
# post-hoc tests with bonf adjustment
summary(glht(maa3, linfct = mcp(trt2 = "Tukey")), test = adjusted("bonf"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = log(amp_acc) ~ trt2 + trialNum + IT_imputed + 
##     (1 + trialNum | colNum), data = sl, REML = TRUE)
## 
## Linear Hypotheses:
##                  Estimate Std. Error z value Pr(>|z|)    
## high - full == 0  0.19453    0.03445   5.647 4.89e-08 ***
## low - full == 0   0.06776    0.02933   2.311   0.0626 .  
## low - high == 0  -0.12677    0.04459  -2.843   0.0134 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)

Generate CI’s for amplitude

# set number of bootstrap samples
nbootSims2 = 1000

# don't need to include hive, but will set trialNum to 2
pframe <- data.frame(trt2 = levels(droplevels(sl$trt2)), 
                     IT_imputed = ITmean,  colNum = 99999, trialNum = 2)
pframe$amp_acc <- 0
# exponentiate to get back to original scale
pp <- exp(predict(maa3, newdata = pframe, re.form=NA, type = 'response')) # re.form sets all random effects to 0


### Calculate CI's (using bootstrap, not accounting for random effects)
bb2 <- bootMer(maa3, FUN=function(x) predict(x, pframe, re.form=NA, type = 'response'), nsim = nbootSims2)
print(paste("Number of bootstrap samples", nrow(bb2$t)))
## [1] "Number of bootstrap samples 1000"
bb2_se <-apply(bb2$t,2,function(x) quantile(x, probs = c(0.025, 0.975)))

#exponentiate to get back to orignal scale
pframe$blo<-exp(bb2_se[1,])
pframe$bhi<-exp(bb2_se[2,])
pframe$predMean <- pp
pframe <- pframe[, c('trt2', "blo", "bhi", "predMean")]
pframe
##   trt2      blo      bhi predMean
## 1 full 39.05177 45.01982 41.91285
## 2 high 46.44376 56.04234 50.91299
## 3  low 41.17250 48.98075 44.85125

Make plots for amplitude for paper

#plot 95% confidence intervals
# "Mean and bootstrap CI based on fixed-effects uncertainty ONLY"
pframe$trt3 = mapvalues(pframe$trt2, from = c("full", "high", "low"), 
                        to = c("Full range\n(220 - 450 Hz)", "High range\n(340 - 390 Hz)", "Low range\n(220 - 330 Hz)"))

pframe
##   trt2      blo      bhi predMean                       trt3
## 1 full 39.05177 45.01982 41.91285 Full range\n(220 - 450 Hz)
## 2 high 46.44376 56.04234 50.91299 High range\n(340 - 390 Hz)
## 3  low 41.17250 48.98075 44.85125  Low range\n(220 - 330 Hz)
ga0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
     geom_point()+ 
     labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency range for reward") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
     theme(axis.text.x = element_text(angle = 45, hjust = 1), 
           legend.position = 'none') +
     theme_classic() + 
  annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 56, label=c("a", "b", "c"),
                color="black") 
ga0

ggsave(plot = ga0, filename = file.path(figDir, "SonicationAmpPredsAndCI_unadjusted.pdf"), width = 5, height = 4)


# make plot using bonferroni adjustments

ga0 <- ggplot(pframe, aes(x=trt3, y=predMean))+
     geom_point()+ 
     labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency range for reward") + 
     geom_errorbar(aes(ymin = blo, ymax = bhi), width = 0.1)+
     theme(axis.text.x = element_text(angle = 45, hjust = 1), 
           legend.position = 'none') +
     theme_classic() + 
  annotate(geom="text", x=c(1,2,3), y=c(0, 0, 0) + 56, label=c("a", "b", "a"),
                color="black") 
ga0

ggsave(plot = ga0, filename = file.path(figDir, "SonicationAmpPredsAndCI_BonfAdjusted.pdf"), width = 5, height = 4)

Visualize aggregated data for amplitude

aggdata <- aggregate(sl$amp_acc, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(aggdata)[colnames(aggdata) == "x"] = "amp"

aggdata_sd <- aggregate(sl$amp_acc, by=list(colNum = sl$colNum,trialNum =sl$trialNum, trt = sl$trt), FUN=sd, na.rm=TRUE)
colnames(aggdata_sd)[colnames(aggdata_sd) == "x"] = "amp_sd"

aggdata = merge(aggdata, aggdata_sd)
aggdata = aggdata[order(aggdata$colNum, aggdata$trialNum, decreasing = FALSE), ]
agg_sm = aggdata[aggdata$trialNum <= 2, ]
rownames(agg_sm) = 1:nrow(agg_sm)
agg_sm
##                colNum trialNum    trt      amp    amp_sd
## 1              blue_4        1   full 40.03641 27.005004
## 2              gold_3        1   full 32.58012 13.395404
## 3           goldred_4        1   full 27.49328  8.575236
## 4             green_4        1   full 44.14520 45.654277
## 5             green_4        2 full_2 19.19985 21.543257
## 6              lime_5        1   full 72.53866 33.940021
## 7          limeblue_5        1   full 49.20279 24.837998
## 8          limeblue_5        2 full_2 59.04592 36.444346
## 9          limegold_5        1   full 34.95839 26.390545
## 10        limegreen_5        1   full 52.56061 32.949195
## 11       limeorange_5        1   full 43.28056 22.966005
## 12       limeorange_5        2 full_2 50.75113 25.907474
## 13         limepink_5        1   full 49.57493 19.711454
## 14       limepurple_5        1   full 43.98249 22.702402
## 15       limepurple_5        2    low 39.95722 19.359774
## 16 limepurpleyellow_5        1   full 34.54496 17.151702
## 17          limered_5        1   full 42.27701 20.963674
## 18          limered_5        2    low 40.63711 20.895896
## 19       limesilver_5        1   full 26.24910  7.336471
## 20        limewhite_5        1   full 37.65233 22.862123
## 21        limewhite_5        2 full_2 40.91118 19.936507
## 22       limeyellow_5        1   full 43.51579 23.156058
## 23           orange_3        1   full 32.25152 18.932874
## 24           orange_5        1   full 55.08235 20.197053
## 25       orangeblue_5        1   full 51.97426 23.301810
## 26      orangegreen_5        1   full 58.59866 33.410963
## 27       orangepink_5        1   full 50.98635 28.370137
## 28     orangepurple_5        1   full 34.53312 16.739467
## 29           purple_3        1   full 53.61466 21.935890
## 30          redblue_4        1   full 28.88506 11.312274
## 31          redblue_4        2 full_2 48.61796 25.277613
## 32         redgreen_5        1   full 78.72927 46.273733
## 33         redgreen_5        2 full_2 70.38830 47.080737
## 34          redpink_5        1   full 74.20294 37.536798
## 35          redpink_5        2    low 56.85081 31.951114
## 36        redpurple_5        1   full 58.06343 40.818616
## 37        redpurple_5        2   high 53.13637 26.813039
## 38           silver_5        1   full 46.34517 30.280907
## 39            white_4        1   full 47.76100 21.162920
## 40            white_4        2 full_2 46.45015 22.564532
## 41        whiteblue_5        1   full 66.11499 33.532056
## 42        whiteblue_5        2    low 67.00735 29.330027
## 43        whitegold_5        1   full 43.10847 26.148044
## 44        whitegold_5        2    low 45.20369 27.241996
## 45       whitegreen_4        1   full 38.87762 18.184245
## 46      whiteorange_5        1   full 61.63802 39.835205
## 47      whiteorange_5        2    low 45.85399 25.838373
## 48        whitepink_5        1   full 63.53760 35.030181
## 49        whitepink_5        2   high 80.19250 47.237615
## 50      whitepurple_5        1   full 47.77247 35.753931
## 51         whitered_5        1   full 48.37233 25.860700
## 52         whitered_5        2   high 74.79446 45.764818
## 53      whiteyellow_5        1   full 65.50425 32.808944
## 54      whiteyellow_5        2 full_2 56.00923 33.372943
## 55       yellowblue_5        1   full 52.45834 30.688932
## 56       yellowblue_5        2   high 79.36156 41.580725
## 57      yellowgreen_5        1   full 61.58710 33.835105
## 58      yellowgreen_5        2 full_2 79.18968 27.697874
## 59     yelloworange_5        1   full 73.32616 36.109449
## 60     yelloworange_5        2   high 54.48161 24.329362
## 61       yellowpink_5        1   full 52.74199 24.869065
## 62       yellowpink_5        2    low 57.76271 29.939244
## 63     yellowpurple_5        1   full 62.77375 29.437162
## 64     yellowpurple_5        2    low 68.98565 35.865129
## 65        yellowred_5        1   full 62.77866 39.753892
## 66        yellowred_5        2 full_2 47.97187 28.376661
agg_sm[duplicated(data.frame(agg_sm$colNum, agg_sm$trialNum)), ]
## [1] colNum   trialNum trt      amp      amp_sd  
## <0 rows> (or 0-length row.names)
ggplot(agg_sm, aes(x = trt, y = amp, fill = trialNum > 1)) +
  geom_boxplot(alpha = 0.2) + 
  geom_point(aes(color = trialNum>1)) + 
  geom_line(aes(group = colNum))

diffdf <- sapply(unique(agg_sm$colNum), FUN = function(x){
  tmp = agg_sm[agg_sm$colNum == x, ]
  if(nrow(tmp) <= 1)
    diff = NA
  else
    diff = tmp$amp[tmp$trialNum == 2] - tmp$amp[tmp$trialNum == 1]
  return(diff)
})

trtDF = sapply(unique(agg_sm$colNum), FUN = function(x){
  tmp = agg_sm[agg_sm$colNum == x, ]
  ttrs = paste(tmp$trt[tmp$trialNum == 1], tmp$trt[tmp$trialNum == 2], sep = "_")
  return(ttrs)
})

buzzdiffs = data.frame(trtDF, diffdf)

ggplot(buzzdiffs, aes(x = trtDF, y= diffdf)) + 
  geom_boxplot() + 
  geom_point()+ 
  labs(y = "amplitude difference m/s/s")
## Warning: Removed 20 rows containing non-finite values (stat_boxplot).
## Warning: Removed 20 rows containing missing values (geom_point).

agg2 = aggregate(sl$amp_acc, by=list(colNum = sl$colNum, fullTrt = sl$trt == "full", trt = sl$trt), FUN=mean, na.rm=TRUE)
colnames(agg2)[colnames(agg2) == "x"] = "amp"
agg2$trt = as.character(agg2$trt)

diffdf <- t(as.data.frame(t(sapply(unique(agg2$colNum), FUN = function(x){
  tmp = agg2[agg2$colNum == x, ]
  if(nrow(tmp) <= 1)
    return(NA)
  if (length(unique(tmp$trt)) > 2){
    tmp = tmp[tmp$trt != "full_2", ]
  }
  diff =  tmp$amp[!tmp$fullTrt] - tmp$amp[tmp$fullTrt]
  return(diff)

}))))

length(diffdf)
## [1] 43
trtDF = sapply(unique(agg2$colNum), FUN = function(x){
  tmp = agg2[agg2$colNum == x, ]
  if (length(unique(tmp$trt)) > 2){
    tmp = tmp[tmp$trt != "full_2", ]
  }
  ttrs = paste(tmp$trt[tmp$fullTrt], tmp$trt[!tmp$fullTrt], sep = "_")
  return(ttrs)
})

length(trtDF)
## [1] 43
buzzdiffs = data.frame(trtDF, diffdf)
tapply(buzzdiffs$diffdf, INDEX = buzzdiffs$trtDF, mean)
##       full_ full_full_2   full_high    full_low 
##          NA    2.871714    9.194234    3.249068
ggplot(droplevels(buzzdiffs[buzzdiffs$trtDF != "full_", ]), aes(x = trtDF, y= as.numeric(diffdf))) + 
  geom_boxplot() + 
  geom_point()+ 
  labs(y = "amplitude difference m/s/s")

Plot amplitude vs. frequency for initial trial, accounting for size

ggplot(sl[sl$trialNum == 1, ], aes(x = freq, y = amp_acc)) + 
  geom_point(position = position_jitter(width = 3), alpha = 0.2) + 
  theme_classic() + 
  geom_smooth(aes(group = trt2)) + 
  labs( y = expression ("Sonication amplitude "(m~s^{-2})), x = "Frequency (Hz)")
## `geom_smooth()` using method = 'gam'

centerFreq = scale(sl[sl$trialNum == 1, "freq"])

m1 <- lmer(log(amp_acc)~centerFreq + I(centerFreq^2) + I(centerFreq^3) + IT_imputed + (1|beeCol), data = sl[sl$trialNum == 1, ])
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(amp_acc) ~ centerFreq + I(centerFreq^2) + I(centerFreq^3) +  
##     IT_imputed + (1 | beeCol)
##    Data: sl[sl$trialNum == 1, ]
## 
## REML criterion at convergence: 3526.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5361 -0.6344  0.0693  0.6857  2.4239 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  beeCol   (Intercept) 0.08395  0.2897  
##  Residual             0.32877  0.5734  
## Number of obs: 1971, groups:  beeCol, 42
## 
## Fixed effects:
##                 Estimate Std. Error t value
## (Intercept)      2.00589    0.57016   3.518
## centerFreq       0.37124    0.03004  12.357
## I(centerFreq^2) -0.08590    0.01227  -7.003
## I(centerFreq^3) -0.03835    0.01055  -3.636
## IT_imputed       0.43530    0.13776   3.160
## 
## Correlation of Fixed Effects:
##             (Intr) cntrFr I(F^2) I(F^3)
## centerFreq  -0.047                     
## I(cntrFr^2) -0.093  0.081              
## I(cntrFr^3) -0.018 -0.844  0.065       
## IT_imputed  -0.996  0.043  0.071  0.019
plot(m1)

# predict -- using mean IT span
ppdf <- data.frame(centerFreq = sort(unique(centerFreq)), IT_imputed = ITmean, beeCol = 99999, acc_amp = 0)
# exponentiate to get back to original scale
ppdf$predAmp = exp(predict(m1, newdata = ppdf, type = "response", re.form = NA))



cf_Unscaled = ppdf$centerFreq * attr(centerFreq, 'scaled:scale') + attr(centerFreq, 'scaled:center')
ggplot(ppdf, aes(x = cf_Unscaled, y = predAmp))+ 
  geom_line() + 
  labs(x = "Sonication Frequency (Hz)", y = expression ("Sonication amplitude "(m~s^{-2})) ) + 
  geom_point(data = sl[sl$trialNum == 1, ], 
             aes(x = freq, y = amp_acc), 
             alpha = 0.2, position = position_jitter(width =2), pch =16, size = .5) 

ggsave(filename = file.path(figDir, "FreqVsAmp_1stTrial_rawData.pdf"), width = 4, height = 3)

g22 <- ggplot(ppdf, aes(x = cf_Unscaled, y = predAmp))+ 
  geom_line() + 
  labs(x = "Sonication Frequency (Hz)", y = expression ("Predicted sonication amplitude "(m~s^{-2})) ) 
g22

ggsave(filename = file.path(figDir, "FreqVsAmp_1stTrial.pdf"), width = 4, height = 3)
sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.1
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] effects_4.0-0   carData_3.0-0   plyr_1.8.4      multcomp_1.4-8 
##  [5] TH.data_1.0-8   MASS_7.3-47     survival_2.41-3 mvtnorm_1.0-6  
##  [9] sjPlot_2.4.0    lme4_1.1-14     Matrix_1.2-11   reshape2_1.4.2 
## [13] ggplot2_2.2.1  
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-131       pbkrtest_0.4-7     RColorBrewer_1.1-2
##  [4] rprojroot_1.2      tools_3.4.2        TMB_1.7.11        
##  [7] backports_1.1.1    R6_2.2.2           sjlabelled_1.0.5  
## [10] DT_0.2             lazyeval_0.2.1     mgcv_1.8-20       
## [13] colorspace_1.3-2   nnet_7.3-12        tidyselect_0.2.3  
## [16] mnormt_1.5-5       compiler_3.4.2     sandwich_2.4-0    
## [19] labeling_0.3       scales_0.5.0       lmtest_0.9-35     
## [22] psych_1.7.8        blme_1.0-4         stringr_1.2.0     
## [25] digest_0.6.12      foreign_0.8-69     minqa_1.2.4       
## [28] rmarkdown_1.8      stringdist_0.9.4.6 pkgconfig_2.0.1   
## [31] htmltools_0.3.6    pwr_1.2-1          htmlwidgets_0.9   
## [34] rlang_0.1.4        shiny_1.0.5        bindr_0.1         
## [37] zoo_1.8-0          dplyr_0.7.4        magrittr_1.5      
## [40] modeltools_0.2-21  bayesplot_1.4.0    Rcpp_0.12.13      
## [43] munsell_0.4.3      abind_1.4-5        prediction_0.2.0  
## [46] stringi_1.1.6      yaml_2.1.14        merTools_0.3.0    
## [49] snakecase_0.5.1    grid_3.4.2         parallel_3.4.2    
## [52] sjmisc_2.6.2       forcats_0.2.0      lattice_0.20-35   
## [55] ggeffects_0.2.2    haven_1.1.0        splines_3.4.2     
## [58] sjstats_0.12.0     knitr_1.17         codetools_0.2-15  
## [61] stats4_3.4.2       glue_1.2.0         evaluate_0.10.1   
## [64] modelr_0.1.1       nloptr_1.0.4       httpuv_1.3.5      
## [67] gtable_0.2.0       purrr_0.2.4        tidyr_0.7.2       
## [70] assertthat_0.2.0   mime_0.5           coin_1.2-1        
## [73] xtable_1.8-2       broom_0.4.2        survey_3.32-1     
## [76] coda_0.19-1        tibble_1.3.4       arm_1.9-3         
## [79] glmmTMB_0.1.4      bindrcpp_0.2